Celem tej pracy domowej jest przeprowadzenie eksploracji zadanego zbioru danych. Dotyczy on pożarów lasów w parku przyrodu Montesinho w północno-wschodniej Portugalii. Dane pochodzą z roku 2007 z publikacji pt. A Data Mining Approach to Predict Forest Fires using Meteorological Data, a dotyczą lat 2000-2003. W tejże pracy możemy znaleźć między innymi mapkę parku, którego dotyczą dane wraz z podziałem jego obszaru na segmenty. 
# Import pakietów
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
Wstępne spojrzenie i zapoznanie się ze zbiorem¶
# Wczytanie danych
fires_df = pd.read_csv("https://lovespreadsheet-tutorials.s3.amazonaws.com/APIDatasets/forest_fires_dataset.csv")
fires_df
Na pierwszy rzut oka zbiorek wygląda w porządku - nie widzimy żadnych NaNów, brakujących wartości czy czegoś, co mogłoby nas mocno zdziwić (może oprócz zerowej powierzchni area - później się temu przyjrzymy).
Teraz wykorzystajmy bazowe narzędzia i zbadajmy dokładniej zbiorek.
fires_df.info()
Okazuje się, że rzeczywiście nie mamy żadnych NaNów w całym zbiorze, ma on 517 wierszy i 13 kolumn, spośród których znaczna część (11) ma wartości liczbowe. Pozostałe dwie kolumny też możemy zamienić na liczbowe - dotyczą one dnia tygodnia i miesiąca odnotowanego pożaru.
Część kolumn ma jednak zdecydowanie enigmatyczne nazwy, dlatego należy sprawdzić, co one oznaczają. W tym celu skorzystamy z dołączonego do zbioru pliku z opisem poszczególnych kolumn.
attributes = pd.read_csv("attributes_forest_fires.csv")
pd.set_option('max_colwidth', 800)
attributes
X i Y to koordynaty obszarów, które zilustrowane są na załączonej na górze mapce.
Wciąż nie wiemy wszystkiego, ale po wiedzę, czym jest system FWI i jego składowe możemy zajrzeć do Internetu albo wprost do wspomnianej wyżej pracy lub też do polskojęzycznej pracy o podobnej tematyce - tutaj.
FWI to Fire Weather Index i jest systemem/standardem do oceny zagrożenia pożarowego. Wsród jego składowych wyróżnia się:
- FFMC - Fine Fuel Moisture Code - oznacza wilgotność górnej warstwy dna lasu - ściółki leśnej, wskaźnik zapalności, najbardziej zmienna składowa;
- DMC - Duff Moisture Code - oznacza wilgotność śrendiej warstwy dna lasu, ma większy czas opóźnienia; również czym większa wartość indeksu, tym niższa wilgotność materiału palnego i wyższy stopeiń zagrożenia;
- DC - Drought Code - oznacza wilgotność dolnej warstwy dna lasu, najwolniej zmieniający się składnik;
- ISI - Initial Spread Index - to wskaźnik początkowej, spodziewanej prędkości propagacji (rozprzestrzeniania się) ognia.
To też dobry moment, żeby sprawdzić, dlaczego w ramce danych z pożarami są pożary, które mają powierchnię 0 ha. Z pracy dowiadujemy się, że te wiersze oznaczają pożary, których powierzchnia była mniejsza niż $\frac{1}{100}$ ha, czyli 100 m$^2$.
Teraz mamy już jakąś szerszą wiedzę i możemy zacząć przyglądać się samym danym...
Podstawowe wielkości opisujące dane¶
fires_df.describe()
Już tutaj możemy zauważyć, że w kolumnach FFMC, ISI czy area są outliery - obserwacje krańcowe znacząco różnią się od kwartyli. Ponadto w opisie zbioru znajdujemy informacje, że zalecana jest transformacja kolumny area ze względu na bardzo dużą skośność. Można to zrobić wykorzystując funkcję $y=log(1+x)$.
Przyjrzyjmy się jednak rozkładom wszystkich zmiennych, może się okazać, że jeszcze któreś wymagałyby transformacji w celu wytrenowania lepszego modelu.
Rozkłady zmiennych¶
fires_df.hist(bins = 25, figsize=(18, 12))
plt.show()
Co widzimy?
Rozkład temperatury temp przypomina normalny, czego moglibyśmy się nawet spodziewać.
Wilgotność powietrza RH, prędkość wiatru wind i parametry pochodzące z Fire Weather Indexu: DMC, DC, ISI również wydają się mieć odpowiednie rozkłady.
Wspominaliśmy już o konieczności zlogarytmowania powierzchni area, natomiast można by się zastanowić nad transformacją również parametru FFMC i opadów deszczu rain.
Przyjrzyjmy się jednak bliżej opadom deszczu...
fires_df["rain"].value_counts()
Okazuje się, że aż 509 rekordów jest powiązanych z dniami, gdzie nie padał deszcz, a tylko 8 dni to takie, gdzie występowały opady. Zatem wydaje się, że nie ma sensu transformacja w tym przypadku.
Spójrzmy jak istotnie różne są histogramy zmiennej area przed i po transformacji, o której wspomnieliśmy. Jako że nie będziemy trenować modelu, nie transformuję na stałe tejże zmiennej.
fig, axes = plt.subplots(nrows=1, ncols=2, figsize=(14, 5))
axes[0].hist(fires_df["area"])
axes[0].set_title("Histogram zmiennej `area` przed transformacją ")
axes[1].hist(np.log1p(fires_df["area"]))
axes[1].set_title("Histogram zmiennej `area` po transformacji y = log(1+x)")
plt.show()
Korelacje i zależności¶
Mamy wiele zmiennych liczbowych, popatrzmy więc najpierw na macierz korelacji między nimi.
plt.figure(figsize=(12, 10))
heatmap = sns.heatmap(fires_df.corr(), annot=True, vmin=-1, vmax=1, cmap="BrBG")
plt.show()
Widzimy korelacje pomiędzy składowymi Fire Weather Indexu i temperaturą oraz wilgotnością powietrza.
Oczywiście poszczególne parametry korelują ze sobą, gdyż są powiązane z wilgotnością kolejnych poziomów gleby oraz szybkością rozprzestrzeniania się pożarów. Korelacje z temperaturą i wilgotnością (w przypadku FFMC) wynikają natomiast z wzorów na te parametry. Także intuicyjnie jesteśmy w stanie stwierdzić, że wilgotność powietrza może mieć wpływ na wilgotność ściółki leśnej.
Te obserwacje trzeba mieć na uwadze przy późniejszym wyborze parametrów do trenowania modelu.
Natomiast możemy zauważyć, że sama zmienna mająca być zmienną wyjaśnianą - area nie jest silnie skorelowana z żadną inną zmienną.
Dla pełniejszego obrazu sytuacji w kwestii zależności spójrzmy na wykresy zależności zmiennych ciągłych.
sns.pairplot(fires_df.drop(["X", "Y"], axis=1))
plt.show()
Możemy zauważyć, że korelacje o których wspomnieliśmy już wyżej są widoczne. Ponadto możemy powiedzieć więcej - zależność FFMC od ISI wydaje się mieć charakter logarytmiczny, a DMC od DC liniowy. Podobnie zachowuje się temperatura względem wilgotności powietrza.
Wybiegając w przyszłość i myśląc o modelu, sprawdźmy również jak wyglądają zależności pomiędzy transformowaną zmienną celu - area a pozostałymi zmiennymi.
months_ordered = ["jan", "feb", "mar", "apr", "may", "jun", "jul", "aug", "sep", "oct", "nov", "dec"]
days_ordered = ["mon", "tue", "wed", "thu", "fri", "sat", "sun"]
fires_df['month'] = pd.Categorical(fires_df['month'], ordered=True, categories=months_ordered)
fires_df['day'] = pd.Categorical(fires_df['day'], ordered=True, categories=days_ordered)
fires_df_transformed = fires_df.copy()
fires_df_transformed["log_area"] = np.log1p(fires_df_transformed["area"])
ax1 = sns.pairplot(fires_df_transformed, y_vars="log_area", x_vars=fires_df_transformed.columns.values[:4], diag_kind = None)
ax2 = sns.pairplot(fires_df_transformed, y_vars="log_area", x_vars=fires_df_transformed.columns.values[4:8], diag_kind = None)
ax3 = sns.pairplot(fires_df_transformed, y_vars="log_area", x_vars=fires_df_transformed.columns.values[8:12], diag_kind = None)
plt.show()
Moim zdaniem nie widać prostych zależności, możemy tylko odnotować pewne obserwacje i zilustrować je lepiej - np. rozkład pożarów (ich wielkości i liczby) w ciągu roku.
Wybrane zależności i wykresy¶
Zacznijmy więc od sprawdzenia, kiedy pojawia się najwięcej pożarów.
plt.figure(figsize=(12,6))
sns.countplot(data=fires_df, x="month", color="orangered")
plt.title("Number of fires over months (in 3 years)", fontsize=16)
plt.xlabel("Month", fontsize=16)
plt.ylabel("Number of fires", fontsize=16)
plt.show()
Najwięcej pożarów ma miejsce w sierpniu i wrześniu.
plt.figure(figsize=(12,6))
sns.countplot(data=fires_df, x="day", color="orangered")
plt.title("Number of fires over days (in 3 years)", fontsize=16)
plt.xlabel("Day", fontsize=16)
plt.ylabel("Number of fires", fontsize=16)
plt.show()
Najwięcej pożarów ma miejsce w weekendy, co może mieć związek z czynnikiem ludzkim.
Sprawdźmy, w których fragmentach parku pożary występują najczęściej.
fires_map = fires_df.groupby(["X", "Y"]).size().reset_index(name='count')
fires_map = fires_map.pivot_table(index="Y", columns="X", values="count", fill_value=0)
#Dodajemy wiersz ze współrzędnymi Y=1 i Y=7, gdzie nie występowały pożary (tylko niewielka cześc parku tam sięga)
fires_map.loc[1] = [0] * 9
fires_map.loc[7] = [0] * 9
fires_map.sort_index(inplace=True)
plt.figure(figsize=(12,6))
ax = sns.heatmap(fires_map, cmap="OrRd", annot=True)
ax.xaxis.set_ticks_position("top")
ax.xaxis.set_label_position("top")
plt.title("Number of fires in park segments", fontsize=16)
plt.show()
Spróbujmy teraz nałożyć tę heatmapę na właściwą mapę parku.
import matplotlib.image as mpimg
map_img = mpimg.imread('park-map.png')
plt.figure(figsize=(12,6))
ax = sns.heatmap(fires_map, cmap="OrRd", annot=True, zorder=2, alpha=0.6)
ax.xaxis.set_ticks_position("top")
ax.xaxis.set_label_position("top")
ax.imshow(map_img,
aspect = ax.get_aspect(),
extent = ax.get_xlim() + ax.get_ylim(),
zorder = 1)
plt.show()
Podobną heatmapę możemy przygotować z medianą powierzchni pożarów w danych segmentach parku.
fires_map_2 = fires_df.groupby(["X", "Y"])["area"].mean().reset_index()
fires_map_2 = fires_map_2.pivot_table(index="Y", columns="X", values="area", fill_value=0)
#Dodajemy wiersz ze współrzędnymi Y=1 i Y=7, gdzie nie występowały pożary (tylko niewielka cześc parku tam sięga)
fires_map_2.loc[1] = [0] * 9
fires_map_2.loc[7] = [0] * 9
fires_map_2.sort_index(inplace=True)
plt.figure(figsize=(12,6))
ax = sns.heatmap(fires_map_2, cmap="OrRd", annot=True, zorder=2, alpha=0.6, fmt=".1f")
ax.xaxis.set_ticks_position("top")
ax.xaxis.set_label_position("top")
ax.imshow(map_img,
aspect = ax.get_aspect(),
extent = ax.get_xlim() + ax.get_ylim(),
zorder = 1)
plt.show()
Segment o tak dużej średniej powierzchni pożarów jest segmentem, gdzie miał miejsce tylko jeden pożar, co możemy zauważyć porównując dwie mapki. Z ich porównania wynika też, że obszarem najbardziej zagrożonym wydaje się wschodnia część parku (o ile mapa jest zorientowana na północ).
Automatyczna eksploracja danych¶
Zgodnie z wymaganiami dotyczącymi zadania mieliśmy wykorzystać również jeden z gotowych narzędzi do automatycznej eksploracji danych. W swojej pracy domowej użyję pandas_profiling.
from pandas_profiling import ProfileReport
report = ProfileReport(fires_df, title="Forest Fires Dataset Profiling Report")
report.to_notebook_iframe()
Wynikiem jest automatycznie wygenerowany raport dotyczący zbioru danych. Zawarta jest w nim znaczna część obserwacji, o których wspominaliśmy wyżej w notatniku, chociażby wstępne overview daje nam informacje o kształcie ramki danych, brakujących wartościach i typach kolumn. Dostajemy też ostrzeżenie o znaczącej ilości zerowych wartości w kolumnach rain i area. Dodatkowo dowiadujemy się, że 4 wiersze w ramce są zduplikowane.
Ponadto otrzymujemy informacje dotyczące każdej zmiennej - kolumny, z których część jesteśmy w stanie uzyskać także samodzielnie używając metody describe ramki danych. Natomiast więcej dowiedzieć możemy się wyświetlając szczegóły, gdzie znajdziemy także więcej parametrów statystycznych dotyczących zmiennych, histogram, najpopularniejsze wartości i najbardziej ekstremalne.
Dalej otrzymujemy możliwość wygenerowania wykresów zależności między kolejnymi zmiennymi liczbowymi oraz macierzy korelacji różnego typu. Oprócz tego wyświetlane są również informacje o brakujących danych i wycinek z ramki danych.
Wykorzystane narzędzie daje nam więc wiele możliwości i znacznie usprawia proces wstępnej analizy danych - otrzymaliśmy większość rezultatów, do których musieliśmy dochodzić krok po kroku.
Jednak należy zauważyć również jego ograniczenia i wady. Pierwszą zauważalną jest czas działania - raport generuje się dość długo, nawet dla tak niewielkiej ramki danych. Raport ma też jednoznacznie określoną strukturę, nie jesteśmy w stanie go dostosowywać do swoich potrzeb (taką potrzebą mogłaby być np. transformacja kolumny area). Co więcej, nie jesteśmy w stanie generować bardziej złożonych wykresów zależności czy bardziej złożonych analiz.
Jest to zatem dobre narzędzie, które może zastąpić nas w wielu kwestiach, ale nie można polegać tylko na nim.
#Raport można zapisać do pliku HTML
report.to_file(output_file='pandas_profiling_report.html')